{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1.3 Dirichlet boundary conditions\n",
    "\n",
    "This section shows how to solve the Dirichlet problem \n",
    "$$\n",
    "-\\Delta u  = f \\quad \\text{ in } \\Omega \n",
    "$$\n",
    "\n",
    "with a *nonhomogeneous Dirichlet boundary condition* $u|_{\\Gamma_D} = g$ on a boundary part $\\Gamma_D$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The extension technique\n",
    "\n",
    "We use the standard technique of reducing a problem with essential non-homogeneous boundary conditions to one with homogeneous boundary condition using an  extension. The solution $u$ in $H^1$ satisfies \n",
    "\n",
    "$$\n",
    "u|_{\\Gamma_D} = g$$\n",
    "\n",
    "and \n",
    "\n",
    "$$ \n",
    "\\int_\\Omega \\nabla u \\cdot \\nabla v_0 = \\int_\\Omega f v_0\n",
    "$$\n",
    "\n",
    "for all $v_0$ in $\\in H_{0,D}^1 = \\{ v \\in H^1: v|_{\\Gamma_D} = 0\\}$. Split the solution \n",
    "\n",
    "$$\n",
    "u = u_0 + u_D\n",
    "$$\n",
    "\n",
    "where $u_D$ is an extension of $g$ into $\\Omega$.   Then we only need to find $u_0$ in $H^1_{0,D}$ satisfying the homogeneous Dirichlet problem \n",
    "\n",
    "$$\n",
    "\\int_\\Omega \\nabla u_0 \\cdot \\nabla v_0 = \\int_\\Omega f v_0 - \\int_\\Omega \\nabla u_D \\cdot \\nabla v_0 \n",
    "$$\n",
    "\n",
    "for all $v_0$ in $H_{0,D}^1$.\n",
    "\n",
    "#### Issues to consider\n",
    " \n",
    " * How to define an extension $u_D$ in the finite element space?\n",
    " * How to form and solve the system for $u_0$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finite element spaces with Dirichlet conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('bottom', 'right', 'top', 'left')"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import netgen.gui\n",
    "%gui tk\n",
    "from ngsolve import *\n",
    "from netgen.geom2d import unit_square\n",
    "\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))\n",
    "mesh.GetBoundaries()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `unit_square` has its boundaries marked as `left`, `right`, `top` and `bottom`. Suppose we want non-homogeneous Dirichlet boundary conditions on \n",
    "\n",
    "$$\n",
    "\\Gamma_D = \\Gamma_{left} \\cup \\Gamma_{right}.\n",
    "$$\n",
    "\n",
    "Then, we set the space as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=2, dirichlet=\"left|right\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare this space with the one without the `dirichlet` flag:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(133, 133)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fs2 = H1(mesh, order=2)\n",
    "fes.ndof, fs2.ndof    # total number of unknowns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, the `dirichlet` flag did not change `ndof`. In NGSolve the unknowns are split into two groups: \n",
    " * dirichlet dofs (or constrained dofs), \n",
    " * free dofs.\n",
    " \n",
    "The facility `FreeDofs` gives a `BitArray` such that FreeDofs[dof] is True iff dof is a free degree of freedom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "free dofs of fs2 without \"dirichlet\" flag:\n",
      " 0: 11111111111111111111111111111111111111111111111111\n",
      "50: 11111111111111111111111111111111111111111111111111\n",
      "100: 111111111111111111111111111111111\n",
      "free dofs of fes:\n",
      " 0: 00001111000011110000111111111111111111110101011011\n",
      "50: 11111111110110110111111111111111011011011111111111\n",
      "100: 111111111111111111111111111111111\n"
     ]
    }
   ],
   "source": [
    "print(\"free dofs of fs2 without \\\"dirichlet\\\" flag:\\n\",\n",
    "      fs2.FreeDofs())\n",
    "print(\"free dofs of fes:\\n\", fes.FreeDofs())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The space `fs2` without `dirichlet` flag has only free dofs (no dirichlet dofs).\n",
    "\n",
    "* The other space `fes` has a few dofs that are marked as *not* free. These are the dofs that are located on the boundary regions we marked as `dirichlet`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forms and assembly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In NGSolve, bilinear and linear forms are defined independently of the dirichlet flags. Matrices and vectors are set up with respect to all unknowns (free or constrained) so they may be restricted to any group of unknowns later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = fes.TrialFunction()\n",
    "v = fes.TestFunction()\n",
    "\n",
    "a = BilinearForm(fes, symmetric=True)\n",
    "a += SymbolicBFI(grad(u)*grad(v))\n",
    "a.Assemble()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $A=$ `a.mat` is the matrix just assembled, then we want to solve for \n",
    "\n",
    "$$\n",
    "  A (u_0 + u_D) = f \\quad \\Rightarrow \\quad A u_0 = f - A u_D\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "  \\left( \\begin{array}{cc} A_{FF} & A_{FD} \\\\ A_{DF} & A_{DD} \\end{array} \\right) \\left( \\begin{array}{c} u_{0,F} \\\\ 0 \\end{array} \\right) = \\left( \\begin{array}{c} {f}_F \\\\ {f}_D \\end{array} \\right) - \\left( \\begin{array}{cc} A_{FF} & A_{FD} \\\\ A_{DF} & A_{DD} \\end{array} \\right) \\left( \\begin{array}{c} u_{D,F} \\\\ u_{D,D} \\end{array} \\right)\n",
    "$$\n",
    "\n",
    "where we have block partitioned using free dofs ($F$) and dirichlet dofs ($D$) as if they were numbered consecutively (which is typically not true). The first row gives\n",
    "\n",
    "$$\n",
    "A_{FF} u_{0,F} = f_F - [A u_D]_F.\n",
    "$$\n",
    "\n",
    "Hence, we need to:\n",
    "\n",
    "- Construct $u_D$ from $g$.\n",
    "- Set up the right hand side from $f$ and $u_D$.\n",
    "- Solve a linear system which involves only $A_{FF}$.\n",
    "- Add solution: $u = u_0 + u_D$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extending the boundary values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we are given that \n",
    "$$\n",
    "g = \\sin(y) \\qquad \\text{ on } \\; \\Gamma_D.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sin(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We interpolate $g$ on the boundary of the domain and extend trivially to $\\Omega$ as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "gfu.Set(g, BND)\n",
    "Draw(gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The keyword `BND` tells `Set` that `g` need only be interpolated on those parts of the boundary that are marked `dirichlet`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve for the free dofs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to assemble the right hand side of $A_{FF} u_{0,F} = f_F - [A u_D]_F$, namely\n",
    "$$\n",
    "r = f - A u_D.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = LinearForm(fes)\n",
    "f += SymbolicLFI(1*v)\n",
    "f.Assemble()\n",
    "\n",
    "r = f.vec.CreateVector()\n",
    "r.data = f.vec - a.mat * gfu.vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The implementation of \n",
    "\n",
    "$$\n",
    "u = \n",
    "u_D + \n",
    "\\left( \\begin{array}{cc} A_{FF}^{-1} & 0 \\\\ 0 & 0 \\end{array} \\right) r\n",
    "$$\n",
    "\n",
    "by sparse solvers is achieved by the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu.vec.data += a.mat.Inverse(freedofs=fes.FreeDofs()) * r  \n",
    "Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using BVP\n",
    "\n",
    "NGSolve also provides a `BVP` facility, within which the above steps are performed automatically. You provide $A$, $f$, a grid function `gfu` with the boundary condition, and a preconditioner. Then `BVP` solves the problem with non-homogeneous Dirichlet boundary condition and overwrites `gfu` with the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = Preconditioner(a,\"local\")   #<- Jacobi preconditioner\n",
    "#c = Preconditioner(a,\"direct\") #<- sparse direct solver\n",
    "c.Update()\n",
    "BVP(bf=a,lf=f,gf=gfu,pre=c).Do()\n",
    "Redraw()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
